Practical 1: A control chart problem
part a - summary of waiting times
Graph displays the daily waiting times period that the Data covers. It is difficult to notice any seasonality or cyclic patterns in the extreme waiting times.
As shown on graph and on table, Monday, Tuesday and Wednesday are the days of the week tht present the highests extreme values. We can also notice that even if Friday has the highest average waiting time, its maximum value is the lowest.
| Weekday | Min | Q1 | Avg | Q3 | Max |
|---|---|---|---|---|---|
| Monday | 116 | 297 | 446 | 559 | 999 |
| Tuesday | 90 | 262 | 431 | 550 | 1154 |
| Wednesday | 95 | 293 | 471 | 597 | 1182 |
| Thursday | 81 | 217 | 373 | 487 | 938 |
| Friday | 148 | 342 | 485 | 622 | 918 |
part b - normal approximation
As the one year period corresponds to \[n = 250\], we can deduce the upper control limit, which corresponds to \[p = 1-1/n\] and equals to 0.996. After scaling this value, we can retrieve the associated quantile or the upper control limit.| x |
|---|
| 987 |
The upper control limit corresponds to 987.472.
Graph suggest that the normal distribution is not appropriate to predict extreme values, the extreme observations tend to diverge strongly from their theorical values under the assumption of a normal distribution. This suggests that the tails are heavy. 
Running the Shapiro-Wilk test shows that the hypothesis of normality is rejected at level 6.10510^{-10}
part c - analzying long waiting times
Block Maxima
The Block Maxima approach divides the dataset into blocks with the same number of observations. Then, for each block we select the observation with the highest value. Looking at the dataset, we can divide it by months or on a weekly basis :
Looking at the data aggregation per month and week, we can say that the first seems to be better. However in order to confirm this hypothesis, when modelling we are going to analyse both cases and keep the model that provides the best results. For the models, we are going to fit a GEV functin on the maxima of each block, represented by the red dots in the following plots.
Peaks-over-threshold approach
The Peaks-over-threshold approach defines the extreme values as the observations above a certain threshold u. In order to find this threshold, we can use an mrlplot and the tcplot.
Looking at the plot, we believe that the best threshold is at a value of 800. However, we are going to see if this holds true by analyzing the results in case we use a threshold of 600.
We are going to fit a GPD on the difference between the observations highleted in the red dots and the thresholds. Once, using a threshold of 800 and then, of 600.
part d - modelling the extremes
Block Maxima
To see what kind of distribution to use, we have to estimate $, $, and $for the maxima. First, we will fit a model on the monthly maxima, then on the weekly ones.
Monthly :
#>
#> Call: fgev(x = maxima_monthly$Average.Wait.Seconds)
#> Deviance: 284
#>
#> Estimates
#> loc scale shape
#> 831.743 144.635 -0.202
#>
#> Standard Errors
#> loc scale shape
#> 38.518 30.520 0.273
#>
#> Optimization Information
#> Convergence: successful
#> Function Evaluations: 119
#> Gradient Evaluations: 87The estimates of the GEV show that the latter is a Weibull: indeed, the shape parameter is negative (-0.202).
Looking at the plot, we remark that if on the one hand the probability plot is not well fitted, on the other hand the observations almost perfectly fit the straight line in the qq-plot. Moreover, the Return Level Plot shows a function that is concave: this is a consequence of the negative shape parameter.
We are now going to make the same analysis but this time using the weekly maxima.
#>
#> Call: fgev(x = maxima_weekly$Average.Wait.Seconds)
#> Deviance: 1236
#>
#> Estimates
#> loc scale shape
#> 585.613 170.495 -0.126
#>
#> Standard Errors
#> loc scale shape
#> 19.8480 14.0297 0.0748
#>
#> Optimization Information
#> Convergence: successful
#> Function Evaluations: 66
#> Gradient Evaluations: 38Once again, we have a shape parameter that is negative, meaning that we are still in a Weibull.
Comparing the results between monthly blocks and weekly ones, we can clearly see that the results are better in the latter. Indeed, points are almost perfectly fitted for both the probability and quantile plot, which was not the case when we fitted the GEV to the monthly maxima. Thus, we can conclude that this for the block maxima approach the maxima from weekly blocks should be used.
POT
We are now going to fit a model using the Peak-over-threshold approach. The function we use for this model is fpot, which fits the difference between the observations and the threshold to a GPD. As we saw previously, two thresholds are likely to be interesting for the modeling: the first one is a threhsold u1 of 800, and the second is at 600. Like we did in the previous analysis with the block maxima, we are going to analyse both and select the best model.
The plots above are the results of the analysis using a threhsold of 600. We can see that the output are good, meaning that the model is well fitted to our data. Moreover, we can observe that once again the Return Level Plot has a function that is concave, suggesting a negative shape parameter for the model.
We are now going to do the same analysis using a threshold of 800
Naturally, having a higher threhsold the number of observations is lower. Looking at the results, we can conclude that we should select the threshold 600 since it allows having a better pp-plot and qq-plot.
part e - upper control limit
Block Maxima
Even if our dataset only consider the working days, we should use 365 days to compute the one year return level.
#> [1] 1295After one year, we expect the Average waiting time in seconds to be of 1295.333 seconds. We are going to use this data in order to compute first the level, then the Value-at-Risk which is going to be providing the confidence line.
We expect the confidence line to be at 1303.624.
Practical 2: Demand monitoring, part I
part a - out of stock products
The red line shows when the observed number of sales goes to zero, essentially the extreme of the sales, there are also other instances where it nears to zero however we ignore them for simplicity sake. We run of stock in these cases.
From the graph above, we can also notice that right before going out of stock, the sales are extremely high. When looking at the first block, from day 0 to day 330 which is the first stockout, the sales data are distributes as follows:
| Sales | Day |
|---|---|
| 0 | 330 |
| 0 | 383 |
| 0 | 384 |
| 0 | 696 |
| 0 | 1274 |
| 0 | 1275 |
One can clearly notice that from the first few observations showed, there is a positive trend in the sales, which the vendors were not able to forecast hence the stockout at Day 330.
We can now look at what happened between the first and the second stockout. 
One can notice that following the first stockout there isn’t a particular structure in the sales data, it looks like the sales are random. Moreover, the stockout lasted 2 days. Hence, we can say that the second stockout is probably due to mismanagement following the initial positive trend of the data.
We can now observe what the sales data between day second and thrid time they went out of stock:
Looking at the sales data distribution in the above graph, one can notice that n the one hand the Sales demand until Day 600 seems kind of linear, but on the other hand following this day it looks like there is once again a positive trend.
Going over each block of data is not necessary as the observations above are already conclusive: the vendors go out of stock when there is an increasing demand. The sellers were probably not able to forecast the positive trends and to manage the inventory given the increasing sales, leading them to stockout.
part b - newsvendor model
Newsvendor model.
Knowing the newsvendor model we can easily compute p:
#> [1] 0.909Alternatively, we could have used the Scperf library and the Newsboy function.
#> Q SS ExpC ExpP CV CR FR z
#> 142.04 41.24 50.03 857.21 0.31 0.91 0.99 1.34Hence, the critical fractile will be computed as \[q = F^-1(0.91)\].
part c - critical fractile and Value at Risk
Value-at-Risk (VaR) also known as chance constraint is used in the newsvendor model as a constraint, limiting the probability of particular event (in this case stock-out) happening. In fact, we can say that the critical fractile of the newsvendor model is the equivalent the level of the Value-at-Risk (0.91): we don’t want to go beyond this value to avoid entering the extreme values.
Expected shortfall or conditional or sometimes called conditional value-at-risk , is about ‘how bad can things get?’. More concrectly, one could compute the Value-at-Risk at a level of 0.91, and compute the expected shortfall which indicates the average value of sales when we enter the extremes of the distribution, hence when the sales are extremely high which, as we saw in the beginning of this analysis, is likely to lead to a stockout.
part d - model using the poisson distribution
#> [1] 0.91Using the proposed function, we can fit the Poisson model to our data an compute the estimated μ by using the MLE.
After having fitted a Poisson to our data, we obtain μ = 100.8. This allows us to compute the critical fractile: according to this model, the estimated critical fractile is 114.
part e - peaks-over-thresholds model
We are now going to compute a POT model. First we have to define the threshold.
Thresholds of 150 and 170 seem interesting, we can test for both.
Comparing the pp-plots of the two models, we decide to use the threshold of 150 as it provides better results. Thus, we will use the model with this threshold to compute the Value at Risk.
#> [1] 145part f - binomial backtesting
First we create the test set that we’re going to use as well as setting
Poisson model
#>
#> Exact binomial test
#>
#> data: sales_violations and nrow(sales_test)
#> number of successes = 89, number of trials = 300, p-value
#> <2e-16
#> alternative hypothesis: true probability of success is not equal to 0.091
#> 95 percent confidence interval:
#> 0.25 0.35
#> sample estimates:
#> probability of success
#> 0.3Looking at the binomial test we see that the p-value is reallt low. Thus, we have to reject the null hypothesis, meaning that the observed number of violations is not equal to the expected number of violations.
POT model
#>
#> Exact binomial test
#>
#> data: sales_violations and nrow(sales_test)
#> number of successes = 15, number of trials = 300, p-value =
#> 0.01
#> alternative hypothesis: true probability of success is not equal to 0.091
#> 95 percent confidence interval:
#> 0.028 0.081
#> sample estimates:
#> probability of success
#> 0.05At a level of 5%, we are once again going to reject the null hyopthesis that te observed number of violations is the same as the estimated one. However, we draw a different conclusion at a level of 1%. In this case when looking at the p-value, one can observe that the null hypothesis can’t be rejected.